//////////////////////////////////////////////////
// CYCLIC DECOMPOSITION - MPI 			//
// Date: April 4, 2012 				// 
// Stillinger Weber equations reference: 	//
// http://lammps.sandia.gov/doc/pair_sw.html 	//
// Paper by Stillinger and Weber 		//
//////////////////////////////////////////////////

#include <math.h>
#include <stdio.h> 
#include "mpi.h"

typedef double prec;

int nWaters = 1024; 

const prec SIGMA = 2.3925;
const prec EPSILON = 6.189;
const prec LAMBDA = 23.15;
const prec a = 1.8;
const prec rCut = 4.3065;
const prec GAMMA = 1.2;
const prec nulltheta = 109.47/180*3.1415926;
const prec A = 7.04955627;
const prec B = 0.6022245584;
int nproc, rank; 

// Structure for Configurations 
typedef struct
{
        prec x, y, z;

} Configuration; 

// Distance Computation 
prec computeDistances(Configuration& ri, Configuration& rj) {

        Configuration cp;
        cp.x = ri.x - rj.x;
        cp.y = ri.y - rj.y;
        cp.z = ri.z - rj.z;

        prec dp = cp.x*cp.x + cp.y*cp.y + cp.z*cp.z;
        prec r = sqrt(dp);

        return r;
}


// Angle Computation between the vectors 
prec computeAngles(Configuration& ri, Configuration& rj, Configuration& rk, prec rij, prec rik) {

        Configuration vectij, vectik;

        vectij.x = ri.x - rj.x; vectij.y = ri.y - rj.y; vectij.z = ri.z - rj.z;
        vectik.x = ri.x - rk.x; vectik.y = ri.y - rk.y; vectik.z = ri.z - rk.z;

        prec dotProduct = vectij.x*vectik.x + vectij.y*vectik.y + vectij.z*vectik.z;

       // prec rij = computeDistances(ri, rj);
       // prec rik = computeDistances(ri, rk);

        prec theta = dotProduct/(rij*rik);

        return theta;

}

// Two-body terms of Stillinger-Weber 
prec twobodyInteraction(Configuration& ri, Configuration& rj) {

        int i, j;
        prec P1 = 0.0;

        prec r = computeDistances(ri, rj);

        if (r < rCut) {

                prec f1 = SIGMA/r;
                prec f12 = f1*f1;
                prec f1p = f12*f12;
                prec f2 = SIGMA/(r-a*SIGMA);
                P1 = A*EPSILON*(B*f1p-1)*exp(f2);

        }

        else P1 = 0.0;

        return P1;
}

prec threeBodyInteraction(Configuration& ri, Configuration& rj, Configuration& rk) {

        prec P2 = 0.0;

        prec rij = computeDistances(ri, rj);
        prec rik = computeDistances(ri, rk);
	
       if ((rij < rCut) && (rik < rCut)) {

                prec theta = computeAngles(ri, rj, rk, rij, rik);

                prec fj = (GAMMA*SIGMA)/(rij-a*SIGMA);
                prec fk = (GAMMA*SIGMA)/(rik-a*SIGMA);

                prec d1 = theta - cos(nulltheta);

                P2 = LAMBDA*EPSILON*d1*d1*exp(fj)*exp(fk);

		}	
        
        else P2 = 0.0;
        
        return P2;
}


prec cyclicDecomposition2B(Configuration *rAtom, int nWaters){

        int i, j;
        prec P1 = 0.0;

        for(i = rank%nproc; i<nWaters; i+=nproc){
                for(j= 0; j<nWaters; j++){
                         if (i<j) {
                         P1 += twobodyInteraction(rAtom[i], rAtom[j]);
			}
                     }
               }

        return P1;

}

prec discDecomposition3B(Configuration *rAtom, int nWaters){

        prec P2 = 0;

        int t = nWaters/3;
        int r = nWaters%3;

        if (r != 0) ++t;

	int p = rank%nproc; 

        for(int i=p; i<nWaters; i+=nproc){

           for(int aa=1; aa<t; ++aa) {

                  int j = (i+aa)%nWaters;

                for(int b=3*aa; b<nWaters; ++b) {

                        int k = (i-aa+b)%nWaters;
                        P2 += threeBodyInteraction(rAtom[k], rAtom[i], rAtom[j]); 
                        P2 += threeBodyInteraction(rAtom[i], rAtom[j], rAtom[k]);
                        P2 += threeBodyInteraction(rAtom[j], rAtom[k], rAtom[i]);

                   }

             }
     }
		
	if (r == 0) { 
		
		for (int i = p; i < t; i+=nproc) { 

			int j = i+t; 
			int k = j+t; 
                        P2 += threeBodyInteraction(rAtom[k], rAtom[i], rAtom[j]); 
                        P2 += threeBodyInteraction(rAtom[i], rAtom[j], rAtom[k]);
                        P2 += threeBodyInteraction(rAtom[j], rAtom[k], rAtom[i]);

		}
	}
     
        return P2;
}

void readPositions(char *fname, Configuration *hAtom) {

        FILE *fp = fopen(fname, "r");
        prec dum;

        // Extract Oxygen Positions from Water and save to posCoord
        for(int j=0; j<nWaters; j++){
                fscanf(fp, "%lf %lf %lf\n",&hAtom[j].x, &hAtom[j].y, &hAtom[j].z);
                fscanf(fp, "%lf %lf %lf\n",&dum, &dum, &dum);
                fscanf(fp, "%lf %lf %lf\n",&dum, &dum, &dum);
        }

        // Close the positions files. 
        fclose(fp);
}


int main(int argc, char *argv[]) { 

	prec global_2B=0.0, global_3B=0.0; 

        Configuration *hAtom = (Configuration*) malloc(nWaters*sizeof(Configuration));

	MPI_Datatype AtomType;
	MPI_Init(&argc, &argv);
 
	MPI_Comm_size(MPI_COMM_WORLD, &nproc); 	
	MPI_Comm_rank(MPI_COMM_WORLD, &rank); 

	// Have processor 0 read in the positions
	if (rank == 0) readPositions("CONFIG_1024.npt", hAtom); 
	MPI_Type_contiguous(3, MPI_DOUBLE, &AtomType);
	MPI_Type_commit(&AtomType); 
	 
	
	prec starttime = MPI_Wtime(); 	

	// Broadcast atom positions to all processors	
	MPI_Bcast(hAtom, nWaters, AtomType, 0, MPI_COMM_WORLD); 

	// prec myPotential_2B = cyclicDecomposition2B(hAtom, nWaters);
	prec myPotential_3B = discDecomposition3B(hAtom, nWaters);

	MPI_Reduce(&myPotential_3B, &global_3B, 1, MPI_DOUBLE , MPI_SUM, 0, MPI_COMM_WORLD);

	prec endtime = MPI_Wtime();
	
	if (rank == 0) { 

                printf("DISC: Time   = %g\n", endtime-starttime); 
                printf("DISC: Energy = %g\n", global_3B);
        }
	
	MPI_Finalize();

	return 0; 
}

	
